39tn aiaa/asme:/sae/asee joint Propulsion Conference and Exhibit Hunisville, AL 20-23 July 2003 AIAA-2003-4919 Calculation of Turbine Axial Thrust by Coupled CFD Simulations of the Main Flow Path and Secondary Cavity Flow in an SLI LOX Turbine D. J. Dorney NASA/George C. Marsliall Space Flight Center Space Transportation Directorate MSFC.AL 35812 Bogdan Marcu Ken Tran Scott Sargent The Boeing Company Rocl<etdyne Propulsion and Power 6633 Canoga Avenue Canoga Park, CA 91309-7922 ABSTRACT Each single reusable Space Launch Initiative (SLI) booster rocket is an engine operating at a record vacuum thrust level of over 730,000 Ibf using LOX and LH2. This thrust is more than 10% greater than that of the Delta IV rocket, resulting in relatively large LOX and LH2 turbopumps. Since the SLI rocket employs a staged combustion cycle the level of pressure is very high (thousands of psia). This high pressure creates many engineering challenges, including the balancing of axial, forces on the turbopumps. One of the main parameters in the calculation of the axial force is the cavity pressure upstream of the turbine disk. The flow in this cavity is very complex. The lack of understanding of this flow environment hinders the accurate prediction of axial thrust. In order to narrow down the uncertainty band around the actual turbine axial force, a coupled, unsteady computational methodology has been developed to simulate the interaction between the turbine main flow path and the cavity flow. The CORSAIR solver, an unsteady three- dimensional Navier-Stokes code for turtx)machinery applications, was used to solve for both the main and the secondary flow fields. Turbine axial thrust values are presented in conjunction with the CFD simulation, together with several considerations regarding the turbine instrumentation for axial thrust estimations during test. Copyright ©2003 by the American Institute of Aeronautics and Astronautics, Inc. No copyright is asserted in the United States under Title 17, U.S. Code. The U.S. Government has a royalty-free license to exercise all rights under the copyright claimed herein for Governmental Purposes. All other rights are reserved for the copyright owner. INTRODUCTION In the recent years the harnessing of increased computing power through parallelization of CFD codes has allowed turtx)machinery engineers to extend their analyses to the secondary flow paths of turbomachines. The prediction of unsteady, three-dimensional viscous flow physics using CFD for the main compressor and turbine flow paths has been coupled with the flow through disk and housing cavities in an effort to better model and understand the actual flow in a turbomachine. One of the main motivations for the analysis of the secondary flow paths has been the necessity to understand and quantify the heat transfer and cooling capacities of such flows, as well as the perfonnance impact the designer must pay in order to use such capacities. Early works in addressing cavity flows [1-4] was driven by the need to understand the cooling cavity flows in basic geometries. The early numerical results reported in the literature [such as Refs. 5 and 6, among an exhaustive list of valuable contributions] have been based on 2-D, steady flow models that do not completely account for the complex environment of rotor-stator interaction coupled with the cavity flow field. Recently [7], results of more complex simulatk>ns have been reported. Such simulations are based on the coupling of large and well- validated CFD solvers for simultaneous calculations of the turbine main flow path (e.g., NASA MS- TURBO) and the cavity flow field (e.g., NASA SCISEAL). Special attention is focused on grid onautics and Astronautics AIAA-2(X)3-4919 issues and solving the flow accurately though the seal packages. While a wealth of information pertaining to the flow of cooling streams in the disk cavity geometries typical to air-breathing gas turbine engines has been, and is being, delivered by the contributions of numerous authors, less data is available for the disk cavity geometries typically encountered in the turbopumps equipping liquid propellant rocket engines. Such geometries have the same level of complexity as the gas turbine engines, but a more severe operating environment. Typically, cryogenic turbopump turbines operate in a hot, high-pressure gas environment, separated by inches from the cryogenic fluid flow path of the pump. If the pumped fluid and the turbine working fluid are compatible, small quantities of cold vaporized fluid from the pump secondary flow paths may be allowed to leak into the turbine cavity. The cold net inflow mixes with hot gas in the cavity in a manner that is not currently understood. If the fluids are not compatible (for example, in a liquid oxygen turbopump powered by gaseous hot hydrogen) the turbine cavity and the pump secondary flows must be carefully organized to keep the fluids separated. One way to achieve separation is through an interstage seal that allows each leakage to drain out. The turbine cavity in this case has a net outflow. In the latter example of a design configuration, a significant uncertainty arises for the case of very large turbopumps; namely the magnitude of the axial thrust produced by the turbine disks. The energetic nature of pump operation requires a carefully organized balance between the axial loads produced on the pump and turbine rotating components. The present paper specifically addresses this uncertainty as related to the Space Launch Initiative RS-83 Engine Liquid Oxygen turt)opump via coupled CFD simulations in which both the turbine main flow path and the cavity flow are simultaneously solved in the full three-dimensional, viscous unsteady regime. The cavity net flow, pressure distribution and resulting axial thrust are calculated for three situations, corresponding to a well sealed cavity with minimum mass flow leaking into the drain system, and two cavities with increased levels of wear for the seals. LOX TURBINE AERODYNAMIC DESIGN The oxidizer turbopump (OTP) turbine for the Boeing Rocketdyne RS-83 hydrogen fueled engine is a single stage, subsonic design. The primary design objectives for this turbine are high performance and robust design. To satisfy these conditions the turbine achieves an efficiency of 77% at the design pressure ratio (PR) of 1.22. Turbine design total inlet pressure is 3876 psia. The turbine has a relatively large pitch line diameter of 16.21 in. with a blade height to diameter ratio (L/D) of 8.9% and a throat aspect ratio (L/O) of 4.69. While the pitch diameter of the turbine is large, the operation of the seal package preventing excessive mass flow leakage out of the turbine cavity imposes a low diameter size for the seal. Consequently, the turbine disk offers a large area for the cavity pressure to act upon and produce axial load. The understanding of the cavity flow and the associated pressure distribution becomes a strong requirement for a robust design. NUMERICAL PROCEDURE The governing equations considered in this study are the time-dependent, three-dimensional Reynolds-averaged Navier-Stokes equations. The numerical algorithm used in the computational procedure consists of a time marching, implicit, finite-difference scheme. The procedure is third- order spatially accurate and second-order temporally accurate. The inviscid fluxes are discretized according to the scheme developed by Roe [8]. The viscous fluxes are calculated using standard central differences. An approximate- factorization technique is used to compute the time rate changes in the primary variables. Newton sub- iterations are used at each global time step to increase stability and reduce linearization errors. The equations of motion are extended to turbulent flows using an eddy viscosity formulation. The turbulent viscosity is calculated using the two-layer Baldwin-Lomax algebraic turbulence model [9]. Message Passing Interface (MPI) and OpenMP directives have been implemented into the code to reduce the computation time for large-scale three- dimensional simulations. The use of MPI allows the coupling of different geometric components, such as the turbine main flow path and cavity, in a straightfonward manner. This procedure was previously demonstrated for unsteady flow in a American Institute of Aeronautics and Astronautics AIAA-2003-4919 supersonic turbine nozzles [10]. with straight centerline The computational procedure uses O- and H- type zonal grids to discretize the flow field and facilitate relative motion of the rotating components. The 0-grids are body-fitted to the surfaces of the airfoils and generated using an elliptic equation solution procedure. They are used to properly resolve tfie viscous flow in the blade passages and to easily apply the algebraic turbulence model. The algebraically-generated H-grids are used to discretize the remainder of the flow field, including the turbine cavity. Figure 1 shows an axial cut through the computational grids. The turbine is modeled as a 53 nozzle/53 blade turbine (1/1 computationally) using 111x45x31 (H) and 121x45x31 (O) grids for nozzle, and 111x45x31 (H) and 121x45x31 (O) grids for rotor. The cavity is modeled using a 35x45x1 1 1 (H) grid. The circumferential span of the computational domain is 6.8 ". The computational analysis has been validated on several rocket-engine turbine geometries, including coupled components such as straight centeriine nozzles (e.g., Refs. 10-14). Boundary Conditions The theory of characteristics is used to determine the boundary conditions at the inlet and exit of the computational domain. For subsonic inlet flow the total pressure, total temperature, and the circumferential and radial flow angles are specified as a function of the radius. The upstream running Riemann invariant, is extrapolated from the interior of the computational domain. For subsonic outflow the circumferential and radial flow angles, total pressure, and the total temperature are extrapolated from the interior of the computatignal domain. For the turbine main flow path the totai-to-static pressure ratio is specified at mid-span of the computational exit and the pressure at all other radial locations at the exit is obtained by integrating the equation for radial equilibrium. For the cavity flow, the upper computational domain is coupled with the main flow computational domain at the H-grid area corresponding to the turbine nozzle hub wall. The cavity walls are one stationary and one rotating with so slip boundary conditions enforced on each. Given the relatively large size of the grid, the labyrinth seal is not modeled numerically, as such additional modeling would have rendered the study computationally too expensive. Instead, the degree of wear in the seal is modeled through a pressure drop across the discharge area at the bottom of the cavity, with the largest the pressure drop corresponding to the most extensive wear modeled. The objective of the calculation is to model the cavity flow at progressively larger mass flow rates. Periodicity is enforced along the outer boundaries of the H-grids in the circumferential direction. For viscous simulations, no-slip boundary conditions are enforced along the solid surfaces. It is assumed that the norrrial derivative of the pressure is zero at solid wall surfaces. In addition, a specified (zero) heat flux is held constant in time along the solid surfaces. Figure 2 shows an instantaneous snapshot of the confluence of the two computational domains and a vector flow visualization of the flow pattem as mass flow is captured by the cavity volume. RESULTS The results presented in this report are circumferentially and time averaged unless specified otherwise. Of interest is the radial distribution of pressures p=p(r) for the three different flow regimes considered. Table 1 presents the mass flow values for t>oth the cavity and the turbine flow path for the three regimes denoted as Low Flow (which is the nominal flow as per design), Medium Flow and High Flow. All denoted values refer to the mass flow leaked at the bottom of the cavity. The algebraic sum of mass flows is not conserved from one case to another: the turbine boundary conditions, which are inlet pressure and temperature and inlet total to exit static pressure ratio, together with the coupling of the secondary cavity flow generate slightly different flow mass values entering the turbine nozzle. Thus, there is some overall variation in the turbine mass flow between the different cases. The variation, however, is less than one tenth of one percent. American Institute of Aeronautics and Astronautics AIAA-2003-4919 Table 1 Mass flow values for the three cases Gbm/s) Case Low Flow Med Flow Higli Flow Turbine 257.4 256.7 256.9 Cavity 3.5 4.8 8.1 Figure 3 shows the radial pressure distribution in the cavity for the three cavity mass flows. The average cavity pressure decreases proportional to the pressure drop across the discharge area at the bottom of the cavity and increasing mass flow. The shape of the radial pressure distribution curve also changes slightly. The fast pressure drop at the top of the cavity for all three cases is explained by examining the vector plot of the cavity flow field in Fig. 4. At the upper lip of the cavity the entering flow forms a stagnation area followed by a downward accelerated flow along the rotating disk wall. The "hump" in all three curves at a radius between 5.75 inches and 6 inches is due to the local geometry (locally increased volume) of the cavity. Table 2 lists the axial loads produced on the turbine disk by the three pressure distributions shown in Fig. 3. The values are calculated by it max F = 2;r jp(r)rdr-^l^P,, (1) Rima where p(r) is the radial pressure distribution, Rmm and Rmax are the radii at the top and bottom of the cavity, and Pexu is the turbine discharge pressure Table 2. Axial load on the turbine disk (values in Ibf) Case Low Flow Med Flow High Flow F -17904.1 -19121.7 -21751.9 F% variation ±3.2% ±3.4% ±4.9% The percent variation in axial thrust is due to significant fluctuation in the pressure distribution values for each case. Figures 5a, 5b and 5c show each of the p = p(r) curves from Fig. 3, as well as it's fluctuation range. The authors did not identify a strong correlation between the rotor blade position and the pressure fluctuation. If such a correlation existed, clocking effects for each cavity sector (of which only a 6.8° period section is modeled here) in relation to the rotor- stator ir^tantaneous position may reduce the percentage variations shown |n Table 2. However, no such correlation has been identified, the pressure variations shown in Fig. 5 being a combined effect of variations in pressure at the upper lip of the cavity and compressibility effects. The pressure distribution p = p(r) is never available when testing hardware, the price and complexity of the instrumentation would be prohibitive. A common practice is to model the pressure variation in the cavity based on a swirling flow model and compute a pressure distribution anchored on the pressure value at the lip of the cavity p(rma)^, which is known at design time from the turbine mean flow path analysis: Pir^)-P^r)= — — (2) IglAA where p is the fluid density, £2 is the turbine angular velocity (in radians per second), g is the gravitational constant, and k is what is usually called as the pumping coefficient or slip factor coefficient (3) The pumping coefficient is a way to estimate how much of the disk tangential velocity U is imparted to the fluid velocity, Cu. The true values of the k factor at each radius k=k(r) are shown in Fig. 6. If one were to use these values (smoothed) and compute the radial pressure distribution in the cavity using expression (2), the results would look as shown in Fig. 7. Clearly, at least for the geometry analyzed here, the swirling flow model is not quite accurate. In practice, for turbopump development testing one or several redundant pressure sensors are placed in the cavity at one particular radius. The choice of where to place the sensors depends on the practice of the team, the availability of space to place the sensor itself, or other particular issues. Figure 8 shows the cavity analyzed in this study and four radial locations at which a static pressure sensor may be placed. These are four distinct possibilities. For each of these test configurations, the local pressure p(rprobe) value is used to calculate a value for the k factor using Eqn. (2). Then, also using American Institute of Aeronautics and Astronautics AIAA-2003-4919 Eqn. (2) and the freshly calculated k factor, a radial pressure distribution is calculated, and integrated to obtain the axial disk load. Figure 9 shows the comparison between the radial pressure distributions shown in Fig. 3 and the pressure distributions obtained based on the four static pressure measurement points, for each flow case. With the exception of the top portion, the radial pressure distribution is matched well for the cases when the choice is to place the pressure sensor at mid-to-lower radii rather than higher radii in the cavity. Optimum for the present case is the KP2 location at the 5.75 in radius, which is also the location of the bulkier local volume in the cavity geometry. This observation is valid for the cavity geometry at hand, but a generalization cannot be made without analyzing a large number of other geometries of cavities. Table 3 lists the values for the axial turbine disk loads computed for the pressure distributions shown in Fig. 9 using expression (1). Table 3. Turbine disk axial loads (Ibf). Case Low Flow Med Flow High Flow F(lbf) %err F(lbf) %err F(lbf) %err KP1 -19785 10.5 -20957 9.6 -24538 12.8 KP2 -17611 1.6 -19106 0.1 -21282 2.1 KP3 -17485 2.3 -18467 3.4 -20660 5.0 KP4 -17185 4.0 -18260 4.5., -20495 5.8 Summarv and Conclusions A series of numerical simulations have been periomied for a coupled turbine/cavity geometry. The simulations were performed for three different leakage rates and levels of wear. The predicted results differed from traditional swirl models, but provide valuable insight into the production of the turbine disk axial thrust loads. The application of a CFD analysis to determine sensor placement locations in turbine cavities was demonstrated. 2. Northrop, A., and Owen, J. M., "Heat-Transfer Measurements in Rotating Disk Systems, Part 2: The Roating Cavity with Radial Outflow of Cooling Air." Int. J. Heat and Fluid Flow, Vol.9, pp.27-36, 1988 3. Chew, J.W., and Rogers, R.H., "An Integral Method for the Calculation of Turbulent Forced Convection in a Rotating Cavity with Radial Outflow", Int. J Heat and Fluid Flow, Vol.9 pp.37, 1988. 4. Graber, D.J., Daniels, W,A., and Johnson, B.V., "Disk Pumping Test:, AFWAL-TR-87-2050, 1987. 5. Virr, G.P., Chew, j. W., and Coupland, J., "Application of CFD to Turbine Disk Cavities", J Turtwmachinery, Vot.116, pp.701 -708, 1994. 6. Athavale, M.M., Przekwas, A.J., Hendricks, R.C. and Steinetz, B.M., "Numerical Analysis of Intra Cavity and Power-Stream Flow Interaction in Multiple Gas-Turbine Disk-Cavities", ASME- 95-GT-325, 1995. 7. Athavale, M.M., Steinetz, B.M. and Hendricks, B.M., "Gas Turbine Primary-Secondary Flowpath Interaction: Transiend, Coupled Simulations and Comparison with Experiments", AIAA-2001-3627, 2001. 8. Roe, P. L, "Approximate Riemann Solvers, Parameter Vectors, and Difference Schemes," Journal of Computational Physics, Vol. 43, 1981, pp. 357-372. Baldwin, B. S., and Lomax, H., "Thin Layer Approximation and Algebraic Model for Separated TuriDulent Flow," AIAA Paper 78- 257, Huntsville, AL, January, 1978. 10. Griffin, L W. and Dorney, D. J., "Simulations of the Unsteady Flow Through the Fastrac Supersonic Turbine," ASME Journal of Turbomachinery, Vol. 122, No. 2, April, 2000, pp. 225-233. References 1. Gan, X., Kilic, M., and Owen, J. K, "Flow Between Counterrotating Disks", J. Turbomachinery, Vol. 117, pp.298-305, 1995 11. Domey, D. J., Griffin, L W., and Huber, F., "A Study of the Effects of Tip Clearance in a Supersonic Turtsine," ASME Journal of Turbomachinery, Vol. 122, No. 4, October, 2000, pp. 674-673. American Institute of Aeronautics and Astronautics AIAA-2(X)3-4919 12. Dorney, D. J., Griffin, L W., Huber, F., Sondal<, D. L., "Effects of Endwali Geometry and Staci<ing on Two-Stage Supersonic Turbine Performance," AIAA Joumal of Propulsion and Power, Vol. 18, No. 6, November-December, 2002, pp. 1305-1308. 13. Domey, D. J., Davis, R. L., Sharma, O. P. "Unsteady f^^ultistage Analysis Using a Loosely Coupled Blade Row Approach". Journal of Propulsion and Power, Vol 12, Nr. 2, P274-282, 1996. 14. Sondak, D. L., and Dorney, D. J., "Vortex Shedding in a Turbine Cascade", International Journal of Turbo and Jet Engines, 16, p 107-126, 1999. American Institute of Aeronautics and Astronautics AIAA-2003-4919 Cavity Figure 1 . The Computational grid American Institute of Aeronautics and Astronautics AIAA-2003-49I9 Figure 2. Interface between the turbine main flow path and cavity flow: fluid is being ingested into the cavity. American Institute of Aeronautics and Astronautics AIAA-2003-4919 3370 3380 3390 Pressure - Psi 3400 3410 3420 Figure 3. Radial pressure distribution in tiie cavity. American Institute of Aeronautics and Astronautics AIAA-2003-4919 7 1 ■ • '^.jjfc 6.5 6 -■.■■■■ ■■^S 5.5 5iTT=5 5 ■ p^':r 4.5 .::■■:■■ ■;-■ . 4 -r:r:=: - 3.5 3 L — J 2.5 ■ «■ ■ 0.5 1 2.5- 0.40.60.8 1 1.2 0.20.40.60.8 1 1.2 a) b) c) Figure 4. Cavity flow field for a) Low Flow b) Medium Flow and c) High Flow cases American Institute of Aeronautics and Astronautics AIAA-2003-4919 -T 1 r- AnsM'Pi -^/-V-V- ifiSiih ^® V=^=-?JL-J vSH 'US ilittlk C=>X= pfyvi :WSti. -5-5557 ' w« -X^X ^^m ^^^r;- s;^^ ^Sfflsffi -ScSB^? SfiS -'-"-"-"'"- rV^y^V-- :====ff?2^elli i=Si=; ----=^-J^J= Vi^_--"- -_- --V.---." _I I 1 I L_ 330 ]» 3(011 MX M 3C9 ntsut'N a) b) c) Figure 5. Cavity radial pressure distributions and the associate ranges of variations for a) Low Flow, b) iVIedium Flow and c) High Flow. American Institute of Aeronautics and Astronautics AIAA-2003-4919 (1 C w 5 3 ^ \ 0.05 Low Flow — Med FLOW High Flow I 1 Ml|||llll '" ■ J_ 0.1 0.15 Pumping factor K value 0.2 0.25 Figure 6. Radial distribution of K factor values. K = Cu/U. American Institute of Aeronautics and Astronautics AIAA-2003-4919 7.5 7 6.5 6 5.5 5- Q: 4.5 4-- 3.5 - 3- 2.5 - o; m. — LowRowP(r) — Med FLOW p(r) High Flow p(r) c Low Flow K-P(r) o Med Flow K-P(r) o High Flow K-P(r) _i_ I ! I r ^ ^'^ ^^■^^r6 ■•& <-> _1_ 3330 3340 3350 3360 3370 3380 Pressure -Psi 3390 3400 3410 3420 Figure 7. Comparison between cavity radial pressure distributions as computed from the CFD flow solution and from the swirl flow model. The K factor values from figure 6 have been used, after smoothing. American Institute of Aeronautics and Astronautics KPl KP2 KP3 KP4 AIAA-2003-4919 2.5 mmimmA 'IttiiniunfiiFi 'IIIISSJKKS' iisaisis iiuii iinii ifiu'.'i) AU III II IMIIi'J Wft\nnHii nm: •>fflH«i II nmii \ omm ir im.- L iiijiuil II Jiifi: ' ii Hill iiiiiiii mU'.. 1 iiiii Ml mm>\ I 'mill mrnvK [ iM \\)\r 0.5 1 Figure 8. Pressure sensor probe locations. American Institute of Aeronautics and Astronautics AIAA-2003-4919 c i s ■ u u i 2S U»Ftt»W -1 1 1 1 1 i !— tmrn-H a) b) c) Figure 9. Cavity radial pressure distribution calculated based on the various K factors anchored based on pressure measurement at the four KP points shown in figure 8. American Institute of Aeronautics and Astronautics